Newark_Cycle_Depth_Rank
### Step 1: Data Collection
# Read Lake Cycle excel sheet
library(readxl)
lake_xls_f <- '/Users/andy/Dropbox/TSCreator/TSCreator development/Developers/Andy/Projects/ML-Data Mining/SpectralAnalysis/LakeCycleHomework_AnalysisExcel.XLS'
excel_sheets(lake_xls_f)## [1] "Plot of Input" "Input Data" "Data Worksheet"
## [4] "Spectral Analysis" "Sheet4"lake_xls <- read_excel(lake_xls_f, sheet='Input Data')
# Complete Interval I of Newark Basin
# Interval I = 2768-3370 m; 706 pts
nb_d <- data.frame(lake_xls[8:713,c(16,17,18)])
colnames(nb_d) <- c('Point', 'Depth (m)', 'Depth Rank')
# Depth Vec
dv <- as.numeric(nb_d$`Depth (m)`)
#diff(dv)
max_d <- max(dv)
min_d <- min(dv)
# window of depth
window <- max(dv)-min(dv)
dr <- as.numeric(nb_d$`Depth Rank`)
max_dr <- max(dr)
min_dr <- min(dr)
plot(nb_d$`Depth (m)`, nb_d$`Depth Rank`, xlab='Depth(meters)',
ylab='Depth Rank', type='l', col='darkblue', lwd=2,
main='Input Data for Spectral Analysis'
)
Step 2. De-mean the data
#________________________________________________________________________________#
nb_d.avg <- mean(dr)
nb_d.dm <- as.numeric(nb_d$`Depth Rank`) - nb_d.avg
# plot
plot(nb_d$`Depth (m)`, nb_d$`Depth Rank`, xlab='Depth(meters)',
ylab='', type='l', col='darkblue', lwd=2,
main='Demeaned Input Data(Red Curve)', ylim=c(-4,4)
)
abline(h=nb_d.avg, lty=2, col='darkblue')
lines(nb_d$`Depth (m)`, nb_d.dm, col='red')
abline(h=0, lty=2, col='red')
Step 3. Filter the data
# Filtering - White noise removed
# (3pt moving "Hanning" average)
# 0.5*n + 0.25*(n-1 + n+1)
hanning_ma <- function(df, n) {
x <- 0.5 * df[n] + 0.25*(df[n-1] + df[n+1])
return(x)
}
nb_d.filt <- c()
i <- 1
nr <- length(nb_d.dm)
for (n in 2:(nr-1)) {
h <- hanning_ma(nb_d.dm, n)
nb_d.filt <- c(nb_d.filt, h)
}
# plot
plot(nb_d$`Depth (m)`, nb_d.dm, xlab='Depth(meters)',
ylab='Value', type='l', col='red', lwd=1,
main='Demeaned + Filtered Input Data', ylim=c(-4,4)
)
abline(h=0, lty=2, col='red')
lines((nb_d$`Depth (m)`)[2:(nr-1)], nb_d.filt, col='green')
Step 4. Split-Cosine Tapering
nb_d.tap <- nb_d.filt
i <- 1
nr <- length(nb_d$`Depth (m)`)
n <- length(nb_d.filt)
tap_perc_n <- n * 10/100 # 10% tapering to reduce the effect of the terminations of data at the ends of the window
for (j in 1:tap_perc_n) {
w = j * pi / tap_perc_n
nb_d.tap[j] = nb_d.filt[j] * 0.5 * (1-cos(w))
}
for (j in (n-tap_perc_n):n) {
w = (512-j) * pi / tap_perc_n
nb_d.tap[j] = nb_d.filt[j] * 0.5 * (1-cos(w))
}
# plot
plot(nb_d$`Depth (m)`, nb_d.dm, xlab='Depth(meters)',
ylab='', type='l', col='red', lwd=1,
main='Demeaned + Filtered + Split cosine tapered(10%) Input Data', ylim=c(-8,8)
)
abline(h=0, lty=2, col='red')
lines((nb_d$`Depth (m)`)[2:(nr-1)], nb_d.filt, col='green')
lines((nb_d$`Depth (m)`)[2:(nr-1)], nb_d.tap, col='black')
# plot(nb_d$`Depth (m)`[], nb_d.tap, xlab='Depth(meters)',
# ylab='', type='l', col='red', lwd=1,
# main='Demeaned + Filtered Input Data', ylim=c(-8,8)
# )Step 5. Fourier Spectral Analysis
#________________________________________________________________________________#
plot.frequency.spectrum <- function(X.k,
xlimits=c(0,length(X.k)/2),
xlab='Frequency',
ylab='Amplitude',
density=FALSE,
plot.type='l'
) {
x <- (0:(length(X.k)-1))
plot.data <- cbind(x, Mod(X.k))
# TODO: why this scaling is necessary?
plot.data[2:length(X.k),2] <- 2*plot.data[2:length(X.k),2]
N <- length(x)
freq <- x
amp <- plot.data[,2]
if(density) {
tot <- sum(amp)
amp <- amp/tot
print(amp)
ylim <- c(0, max(amp))
}
else {
tot <- 1
ylim=c(0,max(Mod(plot.data[,2])))
}
plot(freq, amp, lwd=2, main="", type=plot.type,
#xaxp=c(0, N, N/6),
xlim=xlimits, # show only half of the frequency
ylim=ylim,
xlab=xlab, ylab="Amplitude")
period <- (max_dr - min_dr) /freq
plot(period, amp, lwd=2, main="", type=plot.type,
#xaxp=c(0, N, N/6),
#xlim=xlimits, # show only half of the frequency
#ylim=ylim,
xlab=xlab,
ylab="Amplitude")
}
nb_d.fft <- fft(nb_d.tap)
#par(mfrow=c(2,1))
xlab = sprintf('Frequency (cycle per %f meters)', max_d - min_d)
plot.frequency.spectrum(nb_d.fft,
xlim=c(0,200),
xlab=xlab,
ylab='Amplitude'#,
#,plot.type = 'h'
)
With confidence interval
plot.frequency.spectrum(nb_d.fft,
xlim=c(0,200),
xlab=xlab,
ylab='Amplitude'#,
#,plot.type = 'h'
)
n <- length(nb_d.fft)
R.f <- Mod(nb_d.fft)
I.f <- n * (R.f^2)
var.I.f <- var(I.f)
N <- length(R.f)
freq <- 1:N
df <- 1/sqrt(12)
lines(freq, R.f/qchisq(0.975, df), t='l', lty=2, col='green')
par(new=T)
plot(freq, R.f/qchisq(0.025, df), t='l', lty=2, col='red', xlim=c(0, 200)
, xaxt=NULL, yaxt=NULL, axes=FALSE, xlab="", ylab="")
#
# Calculate Amplitude
Amp <- R.f <- Mod(nb_d.fft)
# Calculate Power Spectrum
PowerSpec <- Amp^2
# Take half of it , coz mirrored
PowerSpec <- PowerSpec[1:(length(Amp)/2)]
PowerSpec.han <- c()
i <- 1
nr <- length(PowerSpec)
for (n in 2:(nr-1)) {
h <- hanning_ma(PowerSpec, n)
PowerSpec.han <- c(PowerSpec.han, h)
}
dp <- as.numeric(nb_d$`Depth (m)`)
window_width <- max(dp) - min(dp)
wavenum <- 1:(length(PowerSpec.han))
wavelen <- window_width/wavenum
TotPowerSpec <- sum(PowerSpec)
RelPowerSpec <- PowerSpec.han/TotPowerSpec
# Charts of Wavenumber vs. Smoothed Power
plot(wavenum, PowerSpec.han, t='h',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Power Spectra [Amplitude]^2'
)
plot(wavenum, PowerSpec.han, t='l',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Power Spectra [Amplitude]^2'
)
plot(wavenum, PowerSpec.han, t='h',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Power Spectra [Amplitude]^2',
xlim=c(0,50))
plot(wavenum, PowerSpec.han, t='l',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Power Spectra [Amplitude]^2',
xlim=c(0,50))
# Charts of Wavenumber vs. Relative Power
plot(wavenum, RelPowerSpec, t='h',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2')
# Charts of Wavenumber vs. Relative Power
plot(wavenum, RelPowerSpec, t='l',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2')
plot(wavenum, RelPowerSpec, t='h',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2',
xlim=c(0,120))
# Charts of Wavenumber vs. Relative Power
plot(wavenum, RelPowerSpec, t='l',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2',
xlim=c(0,120))
# Charts of Wavenumber vs. Relative Power
plot(wavenum, RelPowerSpec, t='h',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2',
xlim=c(0,50))
# Charts of Wavenumber vs. Relative Power
plot(wavenum, RelPowerSpec, t='l',
xlab=paste('Wavenumber (cycle/window) window=', window_width, "m"),
ylab='Relative Power Spectra [Amplitude]^2',
xlim=c(0,50))
Wavelength
# )Wavelength vs. Smoothed Power Spectra
plot(wavelen, PowerSpec.han, t='h',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2')
plot(wavelen, PowerSpec.han, t='l',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2')
plot(wavelen, PowerSpec.han, t='h',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2',
xlim = c(0,40))
plot(wavelen, PowerSpec.han, t='h',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2',
xlim = c(0,10), ylim=c(0,2000))
plot(wavelen, PowerSpec.han, t='l',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2',
xlim = c(0,10), ylim=c(0,2000))
plot(wavelen, PowerSpec.han, t='l',
xlab='Wavelength (m)',
ylab='Power Spectra [Amplitude]^2',
xlim = c(0,100), ylim=c(0,10000))
# )Wavelength vs. Relative Power Spectra
plot(wavelen, RelPowerSpec, t='h',
xlab='Wavelength (m)',
ylab='Relative Power Spectra [Amplitude]^2',
xlim=c(0,10))
abline(v=15, lty=2)
d <- data.frame(WaveLength=wavelen, RelativePowerSpectra=RelPowerSpec)
ord <- order(d['RelativePowerSpectra'], decreasing=TRUE)
RelPowerSpec.sorted <- d[ord,]
RelPowerSpec.sorted$`WaveNumber(WindowWidth=601.68m)` <- wavenum[ord]
summary(nb_d)## Point Depth (m) Depth Rank
## Length:706 Length:706 Length:706
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :characterhead(nb_d)## Point Depth (m) Depth Rank
## 1 1 2768.8 0
## 2 2 2769.66 0.0179745
## 3 3 2770.51 0.0469082
## 4 4 2771.36 0
## 5 5 2772.22 0.348179
## 6 6 2773.07 0.0909985tail(nb_d)## Point Depth (m) Depth Rank
## 701 701 3366.21 0.441702
## 702 702 3367.06 0.172589
## 703 703 3367.92 0.002375
## 704 704 3368.77 0
## 705 705 3369.62 0
## 706 706 3370.48 0.270028head(RelPowerSpec.sorted, 30)## WaveLength RelativePowerSpectra WaveNumber(WindowWidth=601.68m)
## 8 75.210000 0.070802595 8
## 35 17.190857 0.051515077 35
## 26 23.141538 0.047056180 26
## 34 17.696471 0.044246146 34
## 9 66.853333 0.042747939 9
## 27 22.284444 0.042074969 27
## 36 16.713333 0.036760480 36
## 7 85.954286 0.035352798 7
## 25 24.067200 0.027436241 25
## 2 300.840000 0.022828097 2
## 33 18.232727 0.022509834 33
## 37 16.261622 0.022111619 37
## 28 21.488571 0.020977379 28
## 1 601.680000 0.018928085 1
## 3 200.560000 0.015032953 3
## 38 15.833684 0.014104991 38
## 4 150.420000 0.011164447 4
## 32 18.802500 0.010246809 32
## 110 5.469818 0.010094840 110
## 131 4.592977 0.008704774 131
## 39 15.427692 0.008236617 39
## 10 60.168000 0.007889947 10
## 130 4.628308 0.007574582 130
## 29 20.747586 0.007294461 29
## 24 25.070000 0.007263327 24
## 125 4.813440 0.006645236 125
## 5 120.336000 0.006486515 5
## 111 5.420541 0.006224745 111
## 40 15.042000 0.006005567 40
## 109 5.520000 0.005961621 109#________________________________________________________________________________#
###
# Method #1: Assuming a peak is a Milankovitch Cycle; computing Sedimentation Rate
#Eccentricity 412, 123, and 95 k.y.
#Obliquity 41 k.y.
#Precession 23, and 19 k.y.
#Long-period Peak a = 300.84m (2)
#Medium-period Peak b = 75.21 m (8)
#Medium-period Peak c = 23.14m (26)
#Short-period Peak d = 17.19m (35)
#Average of adjacent short -period Peaks e & f = 5.4 m (110, 111)
#Short -period Peak g = 4.6m (130, 131)
hm <- c(8, 35, 26, 2, 110, 131)
mil_cycles <- c(412, 123, 95, 41, 23, 19)
eccentricity <- 123
sed_rate <- (RelPowerSpec.sorted$WaveLength[1]/3) / (mil_cycles/1000)
sed_rate## [1] 60.84951 203.82114 263.89474 611.46341 1090.00000 1319.47368# Unit: m/Ma
# a more accurate sedimentation rate can be computed.plot.harmonic <- function(X.k, i, ts, acq.freq, color="red", xlab='time', add=TRUE) {
X.k.h <- rep(0,length(X.k))
X.k.h[i+1] <- X.k[i+1] # i-th harmonic
harmonic.trajectory <- get.trajectory(X.k.h, ts, acq.freq=acq.freq)
if (add == TRUE)
points(ts, harmonic.trajectory, type="l", col=color)
else
plot(ts, harmonic.trajectory, type="l", col=color, xlab=xlab)
return(harmonic.trajectory)
}
# returns the x.n time series for a given time sequence (ts) and
# a vector with the amount of frequencies k in the signal (X.k)
get.trajectory <- function(X.k,ts,acq.freq) {
N <- length(ts)
i <- complex(real = 0, imaginary = 1)
x.n <- rep(0,N) # create vector to keep the trajectory
f <- 0:(length(X.k)-1)
# Equation x_n = \frac{1}{N} sum_{0}^{N-1} (X_k * exp(i * 2 * pi * \frac{k * n}{ N }))
# X.k : amount of frequency k in the signal, each k-th value is a complex number including strength(amplitude) and phase shift
# N : number of samples
# n : current sample
# k : current frequency, between 0 Hz to N-1 Hz
# 1/N : not necessary but it gives the actual sizes of the time spikes
# n/N : percent of time we have gone through
# 2pik: the speed in radians/second
for(n in 0:(N-1)) { # compute each time point x_n based on freqs X.k
t <- n / N
x.n[n+1] <- sum(X.k * exp(i*2*pi*f*t)) / N
}
x.n * acq.freq
}
plot.show <- function(trajectory,
start_time = 0,
time=1,
harmonics=-1,
plot.freq=FALSE,
plot.type='l',
scale=1.5,
xlab='Time',
ylab='') {
acq.freq <- length(trajectory)/time # data acquisition frequency (Hz)
ts <- seq(start_time,(start_time + time)-1/acq.freq,1/acq.freq) # vector of sampling time-points (s)
X.k <- fft(trajectory)
x.n <- get.trajectory(X.k,ts, acq.freq=acq.freq) / acq.freq
if (plot.freq)
plot.frequency.spectrum(X.k, plot.type=plot.type)
max.y <- ceiling(scale*max(Mod(x.n)))
if (harmonics[1]==-1) {
min.y <- floor(min(Mod(x.n)))-1
} else {
min.y <- ceiling(-scale*max(Mod(x.n)))
}
plot(ts,x.n, type="l", lwd=3, ylim=c(min.y,max.y),
xlab=xlab, ylab=ylab)
#abline(h=min.y:max.y,v=0:time,lty=3, lwd=0.25)
#points(ts,trajectory,cex=0.5, col="red") # the data points we know
y <- rep(0, length(X.k))
if (harmonics[1]>-1) {
for(i in 0:length(harmonics)) {
r <- plot.harmonic(X.k, harmonics[i], ts, acq.freq, color=i+1)
y <- y + r
}
}
points(ts, y/scale, t='l', col='orange', lwd=2)
}
amp <- Mod(nb_d.fft)
N <- length(amp)
amp.h <- amp[0:(N/2)]
x <- 0:(N/2-1)
d <- data.frame(x, amp.h)
d.sorted <- d[order(amp.h, decreasing = TRUE),]
head(d.sorted, 20)## x amp.h
## 9 8 135.57578
## 36 35 90.63089
## 27 26 88.92840
## 35 34 85.61107
## 28 27 82.44720
## 37 36 70.09743
## 3 2 68.08022
## 26 25 59.79256
## 2 1 54.02897
## 5 4 51.95751
## 111 110 50.74177
## 39 38 50.56573
## 38 37 49.17176
## 10 9 46.40244
## 29 28 43.74434
## 132 131 42.32149
## 34 33 41.86663
## 33 32 41.14620
## 126 125 39.84090
## 41 40 37.28518Let’s take the 6 major cycles and try to reconstruct the original signal.
hm <- c(8, 35, 26, 2, 110, 131)
plot.show(nb_d.tap, start_time = min_d, time = max_d - min_d, harmonics = hm,
xlab='Depth (meter)', ylab='', scale=0.5)## Warning in xy.coords(x, y, xlabel, ylabel, log): imaginary parts discarded
## in coercion## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercion
## Warning in xy.coords(x, y): imaginary parts discarded in coercionhm_l <- paste(hm, ' cycles')
hm_l <- c(hm_l, 'combined cycle')
legend('topleft', legend=hm_l, col = c(1, 2, 3, 4, 5, 6,'darkblue'),
lty = rep(1,7), cex=0.40)
Output from R multi-taper library
library(multitaper)
Mspec <- spec.mtm(nb_d.tap, deltat = 0.85, jackknife = TRUE,
k=7, nw=4,
log='no', lwd=2,
xlim=c(0,0.25),
plot=FALSE)
x <- Mspec$freq
y <- Mspec$spec
plot(x, y, type='l', lwd=2
,xlim=c(0,0.2)
#,log='y'
)
The frequencies are plotted in a scale between 0 to 0.5 (In our newark data 0 ~ 0 and 0.5 ~ ). The frequency resolution is low now. But we can still see the peaks around the same frequency.
With Confidence interval.
Mspec <- spec.mtm(nb_d.tap, deltat = 0.85, jackknife = TRUE,
k=7, nw=4,
log='no', lwd=2,
xlim=c(0,0.25)
)
plot(x, y, type='l', lwd=2
,xlim=c(0,0.2)
,log='y', main='Amplitudes ploted on logarithmic scale',
xlab=seq(0,200)
)
Output from R spectrum library
Xspec <- spec.pgram(nb_d.tap, log='dB')
ylimits <- c(1/10*min(Xspec$spec), 25*max(Xspec$spec))
plot(Xspec$freq, Xspec$spec, xlim=c(0,0.2), type='l', lwd=2
#, ylim=ylimits
)
df <- Xspec$df
lines(Xspec$freq, Xspec$spec/qchisq(0.975, df), t='l', lty=2, col='green')
lines(Xspec$freq, Xspec$spec/qchisq(0.025, df), t='l', lty=2, col='red')
bw <- Xspec$bandwidth
df <- Xspec$df
chi_lo <- df / qchisq(0.025, df)
chi_up <- df / qchisq(0.975, df)
#lines(Xspec$freq, Xspec$spec / chi_lo, col='green')
#par(new=T)
#plot(Xspec$freq, Xspec$spec * chi_up, xlim=c(0,0.2), col='red', type='l')
#plot(Xspec, xlim=c(0,0.2))
#plot(Xspec$freq, Xspec$spec, xlim=c(0,0.2), type='l', lwd=2)
#a <- qchisq(Xspec$spec * chi_up, df =Xspec$df)
#lines(Xspec$freq, pchisq(Xspec$spec, df=Xspec$df), col='red', xlim=c(0,0.2)) R psd library output
library(psd)## Loaded psd (1.0.1) -- Adaptive multitaper spectrum estimationd <- ts(nb_d.tap, frequency=1)
Pspec <- psdcore(d, ntapaer=2)
plot(Pspec$freq, Pspec$spec, type='l', xlim=c(0,0.2))
spp <- spectral_properties(Pspec[["taper"]], db.ci = TRUE)
psppu <- with(Pspec, create_poly(freq, dB(spec), spp$stderr.chi.upper))
#par(new=T)
plot(psppu$x.x, psppu$y.y, lty=2, type='l',
xlim=c(0, 0.2), col=c('green','red'),
axes=FALSE, xlab='', ylab='')
library(lomb)
data(ibex)
lsp(ibex[2:3],)
lsp(ibex$temp,times=ibex$hours,type='period',ofac=5)
# lynx contains evenly sampled data
lsp(lynx)
lynx.spec <- lsp(lynx,type='period',from=2,to=20,ofac=5)
summary(lynx.spec)## Value
## Time Time
## Data lynx
## n 114
## Type period
## Oversampling 5
## From 2.0035
## To 19.483
## # frequencies 254
## PNmax 31.579
## At period 9.5763
## At frequency 0.10442
## P-value (PNmax) 1.9603e-12lsp(nb_d.dm)
lsp(nb_d.tap)
lsp(nb_d.tap, type='period')
data(ibex)
ibex.spec <- lsp(ibex[,2:3],type='period',from=12,to=36,ofac=10, plot=FALSE)
op <- par(pch=16)
plot.lsp(ibex.spec, main="Periodogram of daily rhythms of Tb in Capra ibex",
cex.lab=1.3,log="", type="b",level=FALSE,xaxt="n")
axis(side=1,at=seq(12,36,by=6))
par(op)
data(lynx)
set.seed(444)
rand.times <- sample(1:length(lynx),30) # select a random vector of sampling times
randlsp(1000,lynx[rand.times],times=rand.times)## Repeats: 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 960 970 980 990 1000
randlsp(x = nb_d.tap)## Repeats: 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 810 820 830 840 850 860 870 880 890 900 910 920 930 940 950 960 970 980 990 1000